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Abstract 

We demonstrate the use of Langevin spin dynamics for studying dynamical properties of an 
archetypical spin glass system. Simulations are performed on CuMn (20% Mn) where we study the 
relaxation that follows a sudden quench of the system to the low temperature phase. The system is 
modeled by a Heisenberg Hamiltonian where the Heisenberg interaction parameters are calculated 
by means of first-principles density functional theory. Simulations are performed by numerically 
solving the Langevin equations of motion for the atomic spins. It is shown that dynamics is 
governed, to a large degree, by the damping parameter in the equations of motion and the system 
size. For large damping and large system sizes we observe the typical aging regime. 
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I. INTRODUCTION 



Spin glasses exhibit exotic dynamical properties such as aging, memory, and rejuvenation, 
which have triggered a lot of research in the past.-"^ The interest in dynamics in these 
systems is mainly motivated by the fact that practically only off-equilibrium properties 
of spin glasses can be observed in experiments. Peculiarities in the dynamical behavior 
are often related to a complex nature of the phase space and the lack of ergodicity in 
spin-glass systems.- In particular, relaxation towards equilibrium is not characterized by a 
single timescale but rather by a broad spectrum of relaxation times leading to a non-trivial 
evolution of measurable quantities, like e.g. magnetization. 

An aging experiment is an elegant way to reveal the multiscale nature of the spin-glass 
dynamics.^ A system is prepared by quenching from a high temperature to a given tempera- 
ture T and perturbing the system in some way, usually by switching on an external magnetic 
field. A measurement of the time evolution of the magnetization is performed after a certain 
time, tw, has passed since the system preparation. For temperatures T below the spin-glass 
transition temperature Tg, relaxation of the magnetization towards the equilibrium value 
shows a strong dependence on the waiting time, t^. A number of phenomenological mod- 
els has been proposed to explain this behavior, among which are the well-known droplet 
model- and a class of hierarchial models.-'^'^i^iii^ However, scaling laws resulting from these 
models does not allow to interpret experimental results unambiguously, since this would 
require the access to asymptotic regimes of relaxation and hence enormously large time 
scales. Moreover, different, sometimes even contradicting, models give rise to same scaling 
laws, discrediting their predictions. Studying realistic models seems therefore to be a more 
promising way to elucidate mechanisms underlying spin-glass dynamics. 

There has been much work on the non-equilibrium behavior of Edward-Anderson (EA) 
spin-glass models.— 'i^'i^'^'i^'i^ii^'i^ These numerical studies have shown that the dynamical 
behavior of the correlation function is very similar to that of the magnetization (related to the 
correlation function via the fluctuation-dissipation theorem in equilibrium) in experiments. 
In particular, the two-stage relaxation depending on the waiting time is reproduced in the 
simulations. 

In this paper, we study the spin-glass dynamics of the Cu8oMn2o alloy by modeling it 
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with the random-site Heisenberg Hamiltonian: 

^ = - ^ JijCiCjUii ■ nij, (1) 

where rrij represent vector magnetic moments and Cj are the occupation numbers of the 
magnetic atoms (q is equal to 1 if a site is occupied by a Mn atom and other- 
wise); exchange parameters Jij are calculated accurately within the density function theory 
(DFT) approach. Usually, dynamics of spin glasses is studied by means of Monte-Carlo 
simulations.— i^'i^'i^'i^'i^ Such an approach misses some details of local relaxations of spins 
in their local fields. The influence of the finite rate of these local relaxations on the dynamics 
in general can be investigated by solving the Langevin dynamics equations directly. As a 
matter of fact, as will be shown below, motion of individual spins can have a rather strong 
impact on the aging behavior of the system. 

In the next section, we introduce definitions used throughout the paper and describe 
briefly the way aging is observed in experiments. In Section IIIH governing equations for the 
atomistic spin dynamics are given along with some details on the implementation. Results 
of the numerical simulations of the CuMn alloy and EA model are presented in Section IIVI 
The focus is on the influence of damping on the spin dynamics. 

II. DEFINITIONS AND THEORY 

In a typical aging experiment, a system is quenched in zero field from high temperature to 
a temperature below Tg. Then the system is aged during a waiting time t^, a small constant 
field h is applied and the time dependence of the magnetization M(t) or susceptibility 
x{t) = M{t)/h is observed. The relaxation process after the quench can be associated with 
three phases: (1) an initial relaxation towards a local quasi-equilibrium state, (2) aging 
dynamics and (3) global equilibration. The latter is achievable only for systems of finite size 

within a time greater than the ergodic time r^rg ~ exp (A^) . Although equilibration is 
of little interest in an experiment, it must be taken into account in numerical simulations 
when one deals with relatively small systems. 

In numerical simulations, it is convenient to work with the autocorrelation function: 

C{t^ + t,t^) = ^Yl ■ "^^(^^ + ^)]av ' (2) 

i 
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where [. . .]av stand for configuration averaging, i.e. averaging over independent runs with 
randomly generated atomic distributions {q}. In the quasi-equihbrium phase, when the 
conditions for the fluctuation-dissipation theorem (FDT) are fulfilled, the autocorrelation is 
related to the response function^: 

Within linear response, the (thermoremanent) magnetization and susceptibility can be 
found in the following way: 

rtni 

M{t^ + t,t^)=h dt'R{t^ + t,t'), (4) 
Jo 

X{uj,t^) = / dt'R{t^,t')e'^^''-'-\ (5) 



Jo 

where h is a small applied magnetic field in a corresponding thermoremanent magnetization 
(TRM) experiment. In quasi-equilibrium, the relaxation of observables does not depend on 
tw and the relation between the magnetization and the autocorrelation function amounts to 

M(t)=Mo-^Ceq(t), (6) 

where Ceq(t) = C(t^ + t,tw) in a tu,-independent regime. In a non-equilibrium situation, 
when the FDT is violated, different stages of the autocorrelation relaxation can be illustrated 
by the relaxation rate function defined as 

S{t^ + t, tw) = j^^^C*(^ii) + 1, til,). (7) 

The relaxation rate peaks at a time ta of the order of the age of the system, i.e. ta ~ tw 

Spin djTiamics of a spin glass is essentially non-equilibrium, and the aging behavior 
in particular makes the evolution of the autocorrelation function be strongly dependent 
on the waiting time, i.e. time-translation invariance is violated. However, under certain 
circumstances or during certain time intervals, a condition of time-translation invariance, 
C(tyj + 1, tw) = C(t), holds and the relaxation of the system is said to proceed in the quasi- 
equilibrium regime. This assertion can be considered as a definition of a quasi-equilibrium 
state. In order for the true equilibrium dynamics to take place, the fluctuation-dissipation 
theorem must hold. 

Having been prepared at a temperature T, the system tends to the equilibrium state. 
This state can be characterized by the space correlation function which for spin-glasses is 



calculated in the following way: 

i 

Above the transition temperature for sufficiently large i?, the spatial correlation is given by 

G{R)^R'~'+MR/0. (9) 

where d is the dimension of the system, rj a critical exponent, ^ the correlation length and 
u{x) a scaling function which decays to zero for R/^ ^ oo. In the macroscopic limit, the 
correlation length ^ is finite above Tg but diverges with ^ ~ e~'^ as Tg is approached from 
above, where e = (T — Tg)/Tg is the reduced temperature and z/ is a critical exponent. 

III. DESCRIPTION OF METHOD 

Our simulations are performed using the ASD (Atomic Spin Dynamics) packag which 
is based on an atomistic approach for spin dynamics. We use a parametrization of the 
interatomic exchange part of the Hamiltonian in the form of Eq. [TJ The effect of temperature 
is modeled by Langevin dynamics. Connection to an external thermal bath is modeled with 
a Gilbert-like damping characterized by a damping parameter a. 

The microscopic equations of motion for the atomic moments, nij, in an effective field, 
Bj, are expressed as follows: 

^ = -7m, X [B, + hi{t)] - 7-mi X (m^ x [B^ + hi{t)]). (10) 
dt m 

In this expression 7 is the electron gyromagnetic ratio and bj(t) is a stochastic magnetic 
field with a Gaussian distribution. The magnitude of that field is related to the damping 
parameter, a, in order for the system to eventually reach thermal equilibrium. 
The effective field, Bj, on a site z, is calculated from 

arrij 

where for we use the classical Heisenberg Hamiltonian defined by Eq. [H 



We use Heuns method (for details see Ref. 1201) with a time step size of 0.01 femtoseconds 



for solving the stochastic differential equations. Most of the calculations are performed up 
to a time t = 70ps. 



IV. LANGEVIN SPIN DYNAMICS OF A HEISENBERG SPIN GLASS 



A. Spin dynamics of CuMn 

Spin dynamics simulations are performed for the CuMn alloy with 20% of magnetic atoms 
(Mn). The system is described by the Heisenberg Hamiltonian with spins (magnetic atoms) 
distributed over the fee lattice. The magnetic exchange (for the Heisenberg Hamiltonian) pa- 
rameters are obtained by means of the screened generalized perturbation method (SGPM)^ 
implemented within the exact muffin-tin orbital (EMTO)^ scheme. The coherent potential 
approximation (CPA) and disordered local moments (DLM) are used to treat the disordered 
CuMn alloy in the paramagnetic state properly. With this model and for 20 % Mn, critical 
slowing down occurs close to the experimental freezing temperature, Tg, which is approxi- 
mately 90 K.^'^ In contrast to the EA model, where bonds between atomic spins are random 
and the atomic sites ordered, for CuMn, the magnetic exchange parameters between the Mn 
atoms are site-independent and the atomic sites for these atoms are randomly distributed 
in the lattice. Simulations are performed on a 32 x 32 x 32 fee system with 20 % of the 
atomic sites occupied by Mn. We simulate the relaxation process following a quench from 
completely random spin orientations to 10 K. Averaging is performed over 10 random alloy 
configurations with fixed interatomic exchange parameters. 

The main concern of the current work is to investigate the infiuence of damping on 
the aging behavior and on spin relaxation in general, motivation being that the damping 
parameter is the only parameter of the model not derived from first principles calculations. 

First, it is worth considering two limiting cases: a = and a = oo. The first case is trivial 
and corresponds to an utterly deterministic evolution with the total energy conserved. We 
can call this a " microcanonical dynamics" for brevity. Since there is no coupling to a heat 
bath, the system will never reach equilibrium in this type of dynamics. In case of the infinite 
damping parameter, on the other hand, relaxation of a spin to equilibrium with respect to 
its local magnetic field occurs instantaneously, and spin dynamics becomes equivalent to the 
dynamics of the heat-bath Monte Carlo method.— For intermediate values of the damping 
parameter, we expect a system to cross over from the initial off-equilibrium dynamics to 
the regime of relaxation towards a (quasi-) equilibrium state. The duration of the crossover 
must be dependent on the damping parameter. 
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FIG. 1: Autocorrelation + t, tm) calculated for CuMn after a quench from completely random 
spin orientations to T=10 K with damping a = 0.01. The autocorrelation is presented (from left 
to right in the figure) for the logarithmically spaced waiting times tu, = 0, 1.56 • 10^, 3.13 • 10^, 
6.25 • 10^ 1.25 • 10^, 2.5 • 10^, 5 • 10^, 1 • 10^, and 2 • 10"^ fs. 
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FIG. 2: Same as Fig. [T]but with damping a = 0.0316. 



In Figs. [HH] we show the spin glass dynamics of CuMn for four different damping pa- 
rameters, a = 0.01,0.0316,0.1,0.316 respectively. We plot the autocorrelation function for 
logarithmically spaced waiting times. Note that the time is given in femtoseconds, but the 
time step size in the simulation is 0.01 fs. In all four cases, the autocorrelation for ty^=0 
illustrates the initial dynamics of the system right after the quench. The behavior at short 
times {t ^ 500 fs) is similar for all values of a and is shown in Fig. [5] in log-linear scale. The 
evolution of the autocorrelation function can be described here as a sum of two exponents 
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FIG. 3: Same as Fig. [T]but with damping a = 0.1. 




FIG. 4: Same as Fig. [T]but with damping a = 0.316. 

followed by a slower-than-exponential decay at larger times. That is for t ^ 200 fs, we have 

C(t,0) (1 - A)e"^ + Ae~^, (12) 

where A can be extracted from the the intersection point of straight lines corresponding to 
the second exponent: A = 0.32 ± 0.02. The bump in Figs. [Hll] on the curves for = 
corresponds to the second term in Eq. [T21 The initial slope of the curves, 1/ri , is independent 
of damping (see Fig. E]) and temperature (data not shown) and depends only on the details 
of the Hamiltonian and initial spin distribution. When the initial distribution is random, as 
is the case in present simulations, the drop of the autocorrelation is dominated by the strong 
precessional motion of the atomic spins in rapidly varying effective exchange fields. As the 
directions of the effective fields initially are completely randomly oriented, the angle between 
the atomic spin and its effective field is on average large, resulting in a large precessional 
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Time, fs 

FIG. 5: Autocorrelation C{t,tw = 0) for four values of the damping parameter: a = 0.01 (circles), 
a = 0.0316 (boxes), a = 0.1 (triangles), a = 0.316 (crosses). The dashed line is a linear fit to 
the points for 10 fs< t < 40 fs; the slope is equal to 1/ti of Eq. [T2j The slope of the solid lines 
corresponds to the damping relaxation rate 1/t2 in Eg. 1121 

torque on the atomic spins. The system gradually relaxes by means of a damping torque 
on each atomic spin, with the energy of the system dropping down from a high value of 
the random spin configuration ("high-temperature" phase) to a value close to the average 
energy for T = lOK. 

The subsequent decay of the autocorrelation is associated with the equilibration of spins 
in their local fields. Clearly, the rate of this relaxation, 1/t2, depends strongly on the value 
of the damping parameter and for this reason we refer to it as " damping relaxation" . As the 
rate 1/t2 diminishes with increasing a, the initial damping relaxation becomes more difficult 
to identify. As it is seen in Fig. [5] for a = 0.316, the crossover from the initial stage to non- 
exponential decay is rather smooth and relaxation due to damping is indistinguishable. 

The rate of the damping relaxation affects the behavior of the autocorrelation function for 
waiting times much larger than the value of T2. From Fig. [T]one can see that for a = 0.01, the 
curves fall on top of each other for waiting times up to = 625 fs. This implies that by this 
moment, the decay of the autocorrelation function has become time-translation invariant. 
In fact, it seems that for this value of the damping parameter, the system never enters the 
aging regime and the initial relaxation phase crosses over directly to relaxation to the global 
equilibrium for > 5 fs. On the other hand, at smaller values of T2, the aging behavior 
recovers (see Figs. [3|^ and the spin dynamics is similar to the case of infinite damping. 
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FIG. 6: The calculated spin-spin correlation function, 

(mj • Hij), for a = 0.01 (upper panel) and a = 0.1 (lower panel). The correlation function 
is plotted for four logarithmically spaced waiting times. For a = 0.1, all four curves fall essentially 
on top of each other. 

This implies that the strength of damping is determined by the the ratio of the timescale T2 
of the damping relaxation and the characteristic time of detuning of local fields due to the 
motion of neighboring spins contributing to these local fields. 

To determine to what extent the system has equilibrated, one can look at the evolution of 
the spin-spin correlation function g{rij) = (rrij ■ nij). In Fig. [Ulwe plot g{rij) as a function of 
the distance \rij\ between the spins for different waiting times of the system. The correlation 
function is plotted both for a = 0.01 (upper panel) and a = 0.1 (lower panel). As expected 
from the autocorrelation, g{Tij) is seen to evolve faster for the larger damping. It means that 
for sufficiently large damping the system reaches the quasi-equilibrium phase fast enough 
for the aging regime to establish. 

The non-equilibrium behavior seen at small waiting times for small damping parameter 
values can be detected on the microscopic level by observing the trajectories of randomly 
selected spins. In Fig. [7]we plot the trajectory of a typical atomic spin evolving during 100 fs 
(corresponding to a short time scale), for a = 0.01 and a = 0.1, and for two different waiting 
times. The upper panel shows a trajectory for a = 0.01 after a waiting time of 1.25 ps. There 
is a large degree of precessional motion of the atomic spin, confirming the conclusions which 
were drawn from the autocorrelation that the system is still in the initial relaxation phase 
at this waiting time. The middle panel shows the same system after a waiting time of 5 ps. 
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FIG. 7: The trajectory for one typical atomic spin at = 1.25 ps for a = 0.01 (upper panel), at 
= 5 ps for a = 0.01 (middle panel) and at = 1.25 ps for a = 0.1 (lower panel). 

showing an atomic spin with a much more stable spin direction. The spin is now either 
in equilibrium or on the verge of entering equilibrium, although spin motion is much more 
pronounced here than in the aging regime for the system with a = 0.1. The lower panel 
shows the trajectory of an atomic spin for a = 0.1 at a waiting time of 1.25 ps. The system 
is here in the aging regime, as seen in Fig. [21 and the atomic spin direction is stable on a 
time scale of 100 fs. 

In the aging regime, the autocorrelation is characterized by an initial reduction of the 
autocorrelation on to a plateau, similar to what is seen for systems in equilibrium. A plateau 
is clearly seen in Figs. [IH for larger waiting times. The position of the plateau depends on 
temperature and is related to the spin-glass order parameter. More precisely. 



in the macroscopic limit, and the Edward- Anderson order parameter, q-^A, is defined (again 
in the macroscopic limit) as 



where {■ ■ ■)t is thermal averaging. Following the plateau, or the quasi-equilibrium phase, 
is the aging phase. The crossover from one phase to another occurs at a time, tg, when a 
sudden drop of the autocorrelation takes place. The time ts can be best identified as the 
maximum of the relaxation rate defined by Eq. [71 

The relaxation rate for a = 0.1 and for a few waiting times is plotted in Fig. [HI The 
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FIG. 8: Relaxation rate, S{t), for the simulation in Fig. O The relaxation rate is plotted for the 
logarithmically spaced waiting times = 0, 1.56 • 10^ 3.13 • 10^ 6.25 • 10^ 1.25 • 10^ 2.5 • 10^, and 
5 • 10'^ fs. The inset shows the relationship between the waiting time and the peak of the relaxation 
rate for the waiting times t^, = 1.56 • 10^, 3.13 • 10^, 6.25 • 10^, 1.25 • 10^, and 2.5 • 10^ fs. For the last 
waiting time the largest peak seems to be a numerical artifact. Instead the second largest peak 
was chosen as input for the inset. 

relaxation rate is calculated by performing a derivative with respect to Int of the autocor- 
relation. Note that, the poorly defined peaks at the end of the observation time (t ~ 10^ fs) 
are artifacts of a smearing scheme, used when calculating the derivative and which breaks 
down close to the edge of the observation interval. In the inset we show the position of the 
peak relaxation rate, or tg, with respect to the waiting time. As one can see, tg is slightly 
larger than t^, which is expected for the aging regime in a spin-glass system.— However, the 
total time window used in the simulations does not allow to infer any definite form of the 
dependence. 

B. Spin dynamics of the EA model for weak damping 

To investigate even further spin dynamics for small a, we have performed simulations 
of the EA model for a = 0.01 and for different lattice sizes. Simulations have been per- 
formed on a cubic lattice of different sizes L x L x L, where -L=4, 8, 16, and with random 
nearest neighbor exchange interactions drawn from a Gaussian distribution with a standard 
deviation of 1 mRy. This is typically the order of the exchange interaction. The freez- 
ing temperature, Tg, is expected to be 25 K for this model (0.16 within the dimensionless 
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model).— 

In Fig. [HI we show the calculated autocorrelation for a simulation of the Edwards- 
Anderson model. The simulated process is a relaxation following a quench from completely 
random spin orientations to 10 K (0.063 within the standard dimensionless model). The au- 
tocorrelation with respect to observation time is plotted for several logarithmically spaced 
waiting times. As in the case of the CuMn simulations, averaging was performed over 10 
different bond realizations, and for each bond configuration, 10 simulations with different 
initial random spin distributions and different random number sequences in the Langevin 
equations have been done. 

There are three sets of curves in Fig. [9] for three different system sizes. As it is seen, 
for the choice of the damping parameter {a = 0.01) and system sizes {L = 4,8,16), the 
aging regime is very short or even non-present in these simulations. A global equilibrium is 
reached very soon after the initial relaxation has been accomplished. It is also worth noting 
that with the same damping parameter, there is a strong similarity between the waiting 
time dependence of the auto correlation function for the 16x16x16 EA model (Fig. [9]) and 
the dilute 20x20x20 CuMn alloy (Fig. [1]). Moreover, comparing the curves corresponding 
to the initial phase {t^ = 0), one can see that this initial phase is independent of the system 
size. 

Typically, a spin glass system enters the aging regime as soon as local equilibrium condi- 
tions are being met. The dynamics proceeds by a rearrangement of the magnetic order on 
a length scale corresponding to a time scale of the order of the age of the system. In this 
particular simulation, a pure aging regime can not be identified as the system enters a global 
equilibrium soon after the initial relaxation. In contrast to equilibrium, within the aging 
regime, the autocorrelation should depend on the waiting time and not on the system size. 
For the largest four waiting times in Fig. [HI we see the autocorrelation characterized by an 
initial reduction on to a plateau followed by a large sudden reduction to zero for different 
observation times depending on the size of the system. 

A global equilibrium is reached for finite systems, as the correlation length reaches the 
size of the system. If there is aging dynamics, it can no longer proceed. In equilibrium, 
thermal fluctuations continue to govern the motion of the atomic spins. Since there is 
no energy associated with a global rotation of the system, the autocorrelation is reduced 
to zero in equilibrium. The autocorrelation is now translationally invariant with respect 



13 







10° lo' 10^ 10^ lo" lo' 

Observation time (fs) 



FIG. 9: Autocorrelation C(t^ + t,tw) calculated for the Edwards- Anderson model after a quench 
from completely random spin orientations to T=10 K with a = 0.01. Different colors signify 
different system sizes: 4x4x4 (red), 8x8x8 (blue), 16 x 16 x 16 (black). For each system 
size the autocorrelation is presented (from left to right in the figure) for the logarithmically spaced 
waiting times = 0, 1.25 • 10^, 2.5 • 10^, 5 • 10^, 1 • 10^ and 2 • lO'' fs. 

to the observation time, t. By comparing simulations on systems of different size we see 
that by increasing the size of the system by a factor 2, the equilibration time is increased 
approximately by a factor 3. This is due to the fact that the correlation length grows 
logarithmically in time. 

V. CONCLUSIONS 

The investigation of spin dynamics based on the realistic spin-glass model has been per- 
formed by solving the Langevin equations of motion. The exchange parameters have been 
extracted from first-principles DFT calculations, while the damping parameter has been 
varied to study the influence of damping on the dynamics. The simulations showed that 
below the spin-freezing temperature, Tg, the system exhibits the aging behavior for large 
enough values of the damping parameter, a. In this case, the dynamics is very similar to 
that obtained from corresponding Monte-Carlo simulations. 

At weak damping, however, the behavior is different and can basically be characterized 
by two regimes for small and large waiting times, respectively. For waiting times, ty^, below 
some certain value, the autocorrelation function does not depend on (i.e. it is time- 
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translation invariant) and hence is the same as for = 0. For waiting times above a certain 
hmit, the autocorrelation is also time-translation invariant but is characterized by much 
slower decay. The time-translation invariance of the autocorrelation suggests that a system 
of finite size reaches equilibration faster at weaker damping. As a result, it becomes clear 
that spin dynamics inside moderately sized domains in spin-glasses can be strongly affected 
by damping. 
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